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• Abstract 
Q-r 

In the recent Letter 1] we considered the approach of nonequilibrium pair plasma towards 
thermal equilibrium state adopting a kinetic treatment and solving numerically the relativistic 
Boltzmann equations. It was shown that plasma in the energy range 0.1-10 MeV first reaches ki- 
netic equilibrium, on a timescale t^ < 10 -14 sec, with detailed balance between binary interactions 
such as Compton, Bhabha and M0ller scattering, and pair production and annihilation. Later the 
electron-positron-photon plasma approaches thermal equilibrium on a timescale ith ^5 10 -12 sec, 
with detailed balance for all direct and inverse reactions. In the present paper we systematically 

o ; 

present details of the computational scheme used in I], as well as generalize our treatment, consid- 
ering proton loading of the pair plasma. When proton loading is large, protons thermalize first by 

•rH . 

proton-proton scattering, and then with the electron-positron-photon plasma by proton-electron 
scattering. In the opposite case of small proton loading proton-electron scattering dominates over 
proton-proton one. Thus in all cases the plasma, even with proton admixture, reaches thermal 
equilibrium configuration on a timescale tth ^ 10~ n sec. We show that it is crucial to account 
for not only binary but also triple direct and inverse interactions between electrons, positrons, 
photons and protons. Several explicit examples are given and the corresponding timescales for 
reaching kinetic and thermal equilibria are determined. 



1 



I. INTRODUCTION 



An electron-positron plasma is of interest in many fields of physics and astrophysics. One 
of the crucial quantities in this analysis is the timescale of the thermalization process. In 
the early universe Q^^UIs] during the lepton era, ultrarelativistic electron-positron pairs 
contribute to the matter contents of the Universe. In gamma-ray bursts (GRBs) electron- 



positron pairs play essential role in the dynamics of expansion 



Alt 



Indications exist 



on the pre sence of the pair plasma also in active galactic nuclei [9[, in the center of our 



Galaxy 



10| 1 , around hypothetical quark stars 



to appear in the fields of ultra intense lasers 
diagnostic tool for high-energy plasma [ljj. 

In many stationary astrophysical sources the pair plasma is t 
dynamic equilibrium. A detailed study of the relevant processes 



11 ] . In the laboratory pair plasma is expected 
121 ]. where particle production may serve as a 



lought to 

m 



3e in thermo- 



la., B 

22| and spectra 



radiatiation mechanisms [20j, possible equilibrium configurations |l6|.[2l[|. 
231 ] in an optically thin pair plasma has been carried out. Particular attention has been given 

25 1, pair production and annihilation [2^], relativistic 



to collisional relaxation process 



bremsstrahlung [271]. [2a] . double Compton scattering [291]. [301]. 

An equilibrium occurs if the sum of all reaction rates vanishes. For instance, electron- 
positron pairs are in equilibrium when the net pair production (annihilation) rate is zero. 
This can be achieved by variety of ways and the corresponding condition can be represented 
as a system of algebraic equations [31]. However, the main assumption made in all the above 
mentioned works is that the plasma is assumed to obey relativistic quantum statistics. The 
atter is shown to be possible, in principle, in the range of temperatures up to 10 MeV 



14| , 25] . Our main task is to prove that independently of a wide set of initial conditions, 



thermal equilibruim forms for the phase space distribution functions are recovered during 
the process of thermalization by two body and three body direct and inverse particle-particle 
collisions. 

At the same time, in some cases mentioned above the pair plasma can be optically 
thick. Although moderately thick plasmas have been considered in the literature 22j, only 
qualitative description [l^] , [21 ] is available for large optical depths. Assumption of thermal 
equilibrium is often adopted for rapidly evolving systems such as GRBs without explicit 
proof , {§| , |^ • Then hydrodynamic approximation is usually applied both for leptons 
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and photons. However, particles may not be in equilibrium initially. Moreover, they may 
not reach an equilibrium in rapidly evolving systems such as the early Universe or transient 
events, when the energy is released on a very short timescale. 

In the literature there is no consensus on this point. Some authors considered thermal 

33| . In fact, 



equilibrium as the initial state prior to expansion while others did not 

the detailed study of the pair plasma equilibrium configurations, performed in [21], cannot 
answer this question, because essentially nonequilibrium processes have to be considered. 

Thus, observations provide motivation for theoretical analysis of physical conditions tak- 
ing place in nonequilibrium optically thick pair plasma. Notice that there is substantial 



difference between the ion-electron plasma on the one hand and electron-positron plasma 



on the other hand. Firstly, the former is collisionless in the wide range of parameters 
while collisions are always essential in the latter. Secondly, when collisions are important 
relevant interactions in the former case are Coulomb scattering of particles which are usu- 
ally described by the classical Rutherford cross-section. In contrast, interactions in the pair 
plasma are described by quantum cross-sections even if the plasma itself can be still treated 
as classical one. 

Our study reported in [lj],!^! in the case of pure pair plasma clarified the issue of initial 
state of the pair plasma in GRBs sources. Our numerical calculations show that the pair 
plasma on a timescale t < 10~ 12 sec reach thermal equilibrium prior to expansion, due to 
intense binary and triple collisions. In this paper we present details about the computational 
scheme adopted in [lj and turn to a more general case, the pair plasma loaded with baryons. 
Occurence of the thermalization process and the corresponding timescales are necessary for 
determining the dynamics of GRBs. Thermalization timescales t < 10~ 12 sec are indeed 
necessary in order to relate the observed properties of GRBs to the nature of the source, see 
e.g. [36]. 



In the next Section we give qualitative description of the pair plasma, introducing some 
relevant parameters. In Section 3 we discuss pure pair plasma. In Section 4 pair plasma 
with proton loading is discussed. In Section 5 we describe the computational scheme used 
in our analysis. In Section 6 we present results of numerical computations. Discussion 
and conclusions follow in the last Section. In Appendix A relevant conservation laws are 
recalled. In Appendix B conditions for kinetic and thermal equilibria are formulated, and the 
scheme for determination of temperatures and chemical potentials out of number and energy 
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densities are given. Binary interactions in the pair plasma such as Compton, M0ller and 
Bhabha scatterings, as well as pair creation and annihilation by two photons are discussed in 
Appendix C. In Appendix D Compton and Coulomb scatterings with protons are considered. 
In Appendix E three-body radiative variants of the reactions listed above are given. Cutoff 
scheme for numerical evaluation of emission and absorption coefficients are presented in 
Appendix F. In Appendix G mass scaling of the matrix elements for Coulomb scattering 
between electrons, positrons and protons is discussed. In Appendix H the definition of 
matrix elements and cross-sections adopted in the paper are given. 



II. QUALITATIVE DESCRIPTION OF THE PAIR PLASMA 

First of all we specify the domain of parameters characterizing the pair plasma considered 
in this paper. It is convenient to use dimensionless parameters usually adopted for this 
purpose. 

We consider mildly relativistic pair plasma, thus the average energy per particle e brackets 
the electron rest mass energy 

0.1<-^<10. (1) 

The lower boundary is required for significant concentrations of pairs, while the upper bound- 
ary is set to avoid substantial production of other particles such as muons and neutrinos. 

We define the plasma parameter g = (n_c/ 3 ) -1 , where d = \j = ^y/^- ^ s ^ ne Debye 
length, ks is Boltzmann's constant, e, n_ and T_ are the electron charge, number density 
and temperature respectively, c the is speed of light, 6- = ksT^/ '(mc 2 ) is dimensionless 
temperature, u = ^47re 2 n_/m is the plasma frequency and m is the electron mass. To 
ensure applicability of kinetic approach it is necessary that the plasma parameter is small, 
q <C 1. This condition means that kinetic energy of particles dominates their potential energy 
due to mutual interaction. For the pair plasma considered in this paper this condition is 
satisfied. 

Further, the classicality parameter, defined as x= e 2 /(hv r ) = a//3 r , where h is Planck's 
constant, a = e 2 /(hc) is the fine structure constant, v r = (3 r c is mean relative velocity of 
particles, see ( 1F12I) in Appendix. The condition x ^> 1 means that particles collisions can 
be considered classically, while for x <C 1 quantum description is required. In our case both 
for pairs and protons quantum cross-sections are used since x < 1. 
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The strength of screening of the Coulomb interactions is characterized by the Coulomb 
logarithm A = A4dv r /h, where M. is the reduced mass. For electron-electron or electron- 
positron scattering the reduced mass is just m/2, while for electron-proton or positron-proton 
scattering the reduced mass is just the proton mass M. ~ M; for proton-proton scattering 
A4 ~ M/2. Coulomb logarithm varies with mean particle velocity and Debye length, and it 
cannot be set a constant as is usually done in most of studies of the pair plasma. 

Finally, we consider pair plasma with linear dimensions R exceeding the mean free path 
of photons / = (n_cr) _1 , where a is the corresponding total cross-section. Thus the optical 
depth r = naR ^> 1 is large, and interactions between photons and other particles have to 
be taken in due account. We discuss these interaction in the next Section. 

Note that natural parameters for perturbative expansion in the problem under consider- 
ation are the fine structure constant a and the electron-proton mass ratio m/M. 

III. PURE PAIR PLASMA 

For simplicity we first consider pure pair plasma composed of electrons e~ , positrons 
e + , and photons 7. We will turn to a more general case, including protons p in the next 
Section. We assume that pairs or photons appear by some physical process in the region 
with a size R and on a timescale t < R/c. We further assume that distribution functions of 
particles depend neither on spatial coordinates nor on the direction of momenta. We then 
have fi = fi(e,t), namely we consider isotropic distributions functions in momentum space 
for a spatially uniform and isotropic plasma. 

To make sure that classical kinetic description is adequate we estimate the dimensionless 
degeneracy temperature 



and compare it with the estimated temperature in thermal equilibrium. With our initial 
conditions ([T]) the degeneracy temperature is always smaller than the temperature in thermal 
equilibrium and therefore we can safely apply the classical kinetic approach. Besides, since 
we deal with ideal plasma with the plasma parameter g ~ 10~ 3 it is enough to consider only 
one-particle distribution functions. These considerations justify our computational approach 
based on classical relativistic Boltzmann equation. At the same time the right hand side of 
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Boltzmann equations contains collisional integrals as functions of quantum matrix elements, 
as discussed below and in Appendices C-E. 



Relativistic Boltzmann equations 37| , 38j in spherically symmetric case for which the 



original code is designed 



391 ] are 



ldfi n ( dfi ^ TT dfi 
c at \ or r ofi J op 

where /i = cos??, $ is the angle between the radius vector r from the origin and the particle 
momentum p, U is a potential due to an external force, = i^/c are particles velocities, 
fi{e,t) are their distribution functions, the index i denotes the type of particle, e is its 
energy, and 77? and xf are the emission and the absorption coefficients for the production 
of a particle of type "2" via the physical process labeled by q. This is a coupled system of 
partial-integro-differential equations. For homogeneous and isotropic distribution functions 
of electrons, positrons and photons ([3]) reduces to 

i| = EW-xf/i), (4) 

which is a coupled system of integro-differential equations. In P| we also explicitly neglected 
the Vlasov term, describing collisionless interaction of particles in the mean field, since energy 
density of fluctuations of the electromagnetic field are many orders of magnitude smaller 



than the energy density of particles 



40]. 



Therefore, the left-hand side of the Boltzmann equation is reduced to partial derivative 
of the distribution function with respect to time. The right-hand side contains collisional 
integrals, representing interactions between electrons, positrons and photons. 

As example of collisional integral consider absorption coefficient for Compton scattering 
which is given by 

X c 7 7 = y*dk'dpdp'W k i ) p/; k ,p/ 7 (k,t)/ ± (p,t), (5) 

where p and k are momenta of electron (positron) and photon respectively, dp = 
de±doe±j3±/c 3 , dk' = de'^e'^do^/c 3 and the transition function Wk',p';k,p is related to the 
transition probability differential rfwk',p';k,p P er umt time as 

Wu', P ';U,pdk'dp' = Vdw k >^- Kp . (6) 



The differential probability c/uv, P ' ; k,p = Wk'^'^pdk'dp' is given by (1C3I) in Appendix C. 
Given the momentum conservation one can perform one integration over dp' in ([5]) as 



J dp'^k + p-k'-p') -> 1, 



(7) 



but it is necessary to take into account the momentum conservation in the next integration 



over dk.', so we have 



J de 7 5(e 7 + e± 



(8) 



/ d« + 4) 



|9(4 + 4)/^|^ + e± -<-^ 



l^ + e'J/^l 
where the Jacobian of the transformation is 

1 



— J re 



l-/^b'-b' ± ' 



(9) 



and bi = pi/p, bj = p</j/, b' ± = (/?±e±b± + e 7 b 7 - e'^) / ((3' ± e' ± ) . 
Finally, for the absorption coefficient we have 



do' dp 



d|M,«|W 



J c8 f 7 (k,t)f±(p,t), 



(10) 



16e±e 7 e' ± 

where the matrix element here is dimensionless. This integral is evaluated numerically as 
described in Appendix. 

For all binary interactions we use exact QED matrix elements which can be found in the 
standard textbooks, e.g. in [4l|.[42|.|43l|. and are given in Appendix [Cl 

In order to account for the charge screening we introduced the minimal scattering angles 
following 44] . see Section [F] in Appendix. This allows to apply the same scheme for the 
computation of emission and absorption coefficients for Coulomb scattering, while many 
treatments in the literature use the Fokker-Planck approximation, e.g. 45]. 

For such a dense plasma collisional integrals in (jlj) should include not only binary inter- 



actions, having order a 2 in Feynmann diagrams, but also triple ones, having order a 3 
As example for triple interactions consider relativistic bremsstrahlung 



4l| 



ei + e 2 ^ e[ + e' 2 + 7'. 



(11) 
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For the time derivative, for instance, of the distribution function f 2 in the direct and in the 
inverse reactions (fTTj) one has 

h = J dpidpidpadk' [W p > iiP ' 2 y. >Pl>P2 f[f2f' k - 

-W pi>P2;p > liP > 2>k ff 1 f 2 ] = J rfpidp^p^k'-^-pX (12) 



SW(P f -Pi)\M fi \ 
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2 5 eie 2 e' ie ^' 



where 



dp 1 dp 2 W p ' i>p > 2 y. puP2 = V 2 dw 1 , 
dp[dp' 2 dk'W Pl>pr>p > vP2> v = Vdw 2 , 

and dwi and dw 2 are given by (1H3j) for the inverse and direct process (TTTj) respectively. The 
matrix element here has dimensions of the length squared, see Section [H] in Appendix. 
In the case of the distribution functions (|T5l) . see below, we have multipliers proportional 

to 

Ft = exp ^, (13) 

called fugacities, in front of the integrals. The calculation of emission and absorption coeffi- 
cients is then reduced to the well known thermal equilibrium case [31]]. In fact, since reaction 
rates of triple interactions are a times smaller than binary reaction rates, we expect that 
binary reactions come to detailed balance first. Only when binary reactions are all balanced, 
triple interactions become important. In addition, when binary reactions come into balance, 
distribution functions already acquire the form (fl5l) . Although there is no principle diffi- 
culty in computations using exact matrix elements for triple reactions as well, our simplified 
scheme allows for much faster numerical computation. The corresponding reaction rates for 
triple interactions are given is Section [E] in Appendix. 

We consider all possible binary and triple interactions between electrons, positrons and 
photons as summarized in table [B 

Each of the above mentioned reactions is characterized by the corresponding timescale 
and optical depth. For Compton scattering of an electron, for instance, we have 

t cs = — - — , r cs = cr T n±R, (14) 



S 



Binary interactions 


Radiative and 
pair producing variants 


M0ller and Bhabha 

6?i Co y Ci Co 
1 Z 1 Z 

e ± e T _^ e ±' e T' 


Bremsstrahlung 

1 Z 1 Z ' 


Single Compton 

e ± 7 — >e ± 7' 


Double Compton 

e 7<->e '/'y" 


Pair production 
and annihilation 


Radiative pair production 
and 3-photon annihilation 

77'<— >e =1: e =F 7" 
e =1: e = ' = <-^77 / 7" 







TABLE I: Microphysical processes in the pair plasma. 



where ot = ^f a2 (^) 2 is the Thomson cross-section. There are two timescales in our problem 
that characterize the condition of detailed balance between direct and inverse reactions, t cs 
for binary and a _1 t cs for triple interactions respectively. 

We choose arbitrary initial distribution functions and find a common development. At a 
certain time the distribution functions always have evolved in a functional form on the 
entire energy range, and depend only on two parameters. We find in fact for the distribution 
functions the expressions 



Me) 



exp 



(2vr/i) 3 V Oi 
with chemical potential Vt = and temperature 0; = where e = 
the particle. Such a configuration corresponds to a kinetic equilibrium 



(15) 

is the energy of 



45) 



461 ] in which 



particles acquire a common temperature and nonzero chemical potentials. At the same 
time we found that triple interactions become essential for t > t^, after the establishment of 
kinetic equilibrium. In strict mathematical sense the sufficient condition for reaching thermal 
equilibrium is when all direct reactions are exactly balanced with their inverse. Therefore, 
in principle, not only triple, but also four-particle, five-particle and so on reaction have to 
be accounted for in equation (j3j). The timescale for reaching thermal equilibrium will be 
then determined by the slowest reaction which is not balanced with its inverse. We stress, 



however, that the necessary condition is the detailed balance at least in triple interactions, 
since binary reactions do not change chemical potentials. 



Notice that a method similar to ours was applied in [45j in order to compute spectra of 
particles in kinetic equilibrium. However, although the approach was similar, the computa- 
tion was never carried out in order to actually observe the reaching of thermal equilibrium. 

Finally, it is worth mentioning the physical meaning of the chemical potential i\ in kinetic 
equilibrium entering the formula ffl5l) . In the case of pure pair plasma a non-zero chemical 
potential represents deviation from the thermal equilibrium through the relation 

u k = 9\n(n k /n th ), (16) 

where n t h are concentrations of particles in thermal equilibrium. 

IV. PROTON LOADING 

So far we dealt with leptons, having the same mass but opposite charges. In that case the 
condition of electric neutrality is identically fulfilled. We described electrons and positrons 
with the same distribution function. Situation becomes more complicated when admixture 
of protons is allowed. Since charge neutrality 

n_ = n + + n p (17) 

is required, the number of electrons is not equal to the number of protons. In such 
new dimensionless parameter, the baryonic loading B, can be introduced as 

NMc 2 n p Mc 2 
B = —^ = —> (18) 

where N and n p are the number and the concentration of protons, £ and p r = p 7 +p + +p_ are 
radiative energy and energy density respectively. Since in relativistic plasma electrons and 
positrons move with almost the speed of light, both photons and pairs in thermal equilibrium 
behave as relativistic fluid with equation of state p r ~ p r /3. At the same time, protons are 
relatively particles in the energy range (pQ), with negligible pressure and dust-like equation of 
state p ~ 0. In this way by introducing parameter B we distinguish a radiation-dominated 
(B < 1) from a matter-dominated (B > 1) plasma. For electrically neutral plasmas there 
exists an upper limit on the parameter B defined by (|18p . which is B <M/m. 
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Binary interactions 


Radiative and 
pair producing variants 


Coulomb scattering 


Bremsstrahlung 


P1P2 — > v'iP2 

± ' ±/ 
pe — > p e 


P1P2 <-> p'lP' 2 l 

± ' ±/ 
pe npe 7 




pef <-» p'e^'e^e^ 


Single Compton 


Double Compton 


P7 — >p'j' 


pj^p'j'j" 




pj^-f-p'e^e 11 



TABLE II: Microphysical processes in the pair plasma involving protons. 

In the range of energies ([1]) the radiative energy density can be approximated as p r ~ 
n_mc 2 , and then we have for concentrations n p ~ n_B^. If protons and electrons are at 
the same temperature then from the equality of the kinetic energy of a proton €k tP = ^ E 
and the one of an electron ek- ~ mc 2 we have ^ ~ v/^, therefore protons are indeed 
nonrelativistic. 

In presence of protons additional binary reactions consist of Coulomb collisions between 
electrons (positrons) and protons, scattering of protons on protons and Compton scattering 
of protons. Additional triple reactions are radiative variants of these reactions, see Table HT1 
and Appendix ID1 

Protons can be thermalized in two ways: either in a two-step process first between them- 
selves and then by electron/positron-proton collisions, or just by the latter mechanism. The 
rate of proton-proton collisions is a factor \p%^ L ~ B (jj) 3 ^ 2 smaller than the rate of 
electron-electron collisions, see (ID15I) . The rate of proton-electron/positron collisions is a 
factor j-js ~ H smaller than the one of electron-electron collisions, see (IDllj) . Therefore, for 
B > */|f proton-proton collisions are faster, while for B < proton-electron/positron 
ones predominate. 
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V. THE DISCRETIZATION PROCEDURE AND THE COMPUTATIONAL 
SCHEME 



In order to solve equations (jl]) we use a finite difference method by introducing a com- 
putational grid in the phase space to represent the distribution functions and to compute 



collisional integrals following 39|. Our goal is to construct the scheme implementing energy, 
baryon number and electric charge conservation laws, see Appendix |A] For this reason we 
prefer to use in the code, instead of distribution functions /j, the spectral energy densities 

EM) = (19) 

& 

where = \J\ — (mid 2 / e^) 2 , in the phase space e^. Then 

Eifi(p,t)drdp = —^-drdei = Eidrdei (20) 

is the energy in the volume of the phase space dvdp. The number density of particles of 
type "i" is given by 

n% = J fidp = J ^-dei, drii = fcdp, (21) 
while the corresponding energy density is 

Pi = / tifidp = / Eidei. 



We can rewrite Boltzmann equations (jlj) in the form 

^T = £(fl-xW' ( 22 ) 
g 

where f]\ = (inefPi/c 3 )^. 

We introduced the computational grid for phase space {ej,/i, cp}, where \i = cos^, $ and 
are angles between radius vector r and the particle momentum p. The zone boundaries 
are e i)U , Tl/2 , /i fc=F i/ 2 , 0* T i/ 2 for 1 < u < cu max , 1 < k < /c max , 1 < I < l max . The length of the 
i-th interval is Aej iW = — £i,u~i/2- On the finite grid the functions (Tl9|) become 

= [ deEfe). (23) 



Now we can replace the collisional integrals in (I22p by the corresponding sums. 
After this procedure we get the set of ordinary differential equations (ODE's), instead 
of the system of partial differential equations for the quantities E iiUJ to be solved. There 
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are several characteristic times for different processes in the problem, and therefore our 
system of differential equations is stiff. Under these conditions eigenvalues of Jacobi matrix 
differs significantly, and the real parts of eigenvalues are negative. We use Gear's method 



4=71 ] to integrate ODE's numerically. This high-order implicit method was developed for the 



solution of stiff ODE's. 

In our method exact energy conservation law is satisfied. For binary interactions the 
particles number conservation law is satisfied as we adopt interpolation of grid functions 
E iu) inside the energy intervals. 



VI. NUMERICAL RESULTS 



In what follows we consider in details three specific cases. In the first two cases our grid 
consists of 60 energy intervals and 16 x 32 intervals for two angles d and 4> characterizing 
the direction of the particle momentum. In the third case we have 40 energy intervals. 



A. Case I 



We take the following initial conditions: flat initial spectral densities ^(e^) = const, 
total energy density p = 10 24 erg/cm 3 . Plasma is dominated by photons with small amount 
of electron-positron pairs, the ratio between energy densities in photons and in electron- 
positron pairs ( = p±/p 7 = 10 -5 . Baryonic loading parameter B = 10~ 3 , corresponding to 
p p = 2.7 x 10 18 erg/cm 3 . The energy density in each component of plasma changes, as 
can be seen from fig. [TJ keeping constant the total energy density shown by dotted line in 
fig. dj as the energy conservation requires. As early as at 10~ 23 sec the energy starts to 
redistribute between electrons and positrons from the one hand and photons from the other 
hand essentially by the pair-creation process. This leads to equipartition of energies between 
these particles at 3 x 10~ 15 sec. Concentrations of pairs and photons equalize at 10 -14 sec, 
as can be seen from fig. [2l From this moment temperatures and chemical potentials of 
electrons, positrons and photons tend to be equal, see fig. [3] and fig. H] respectively, and it 
corresponds to the approach to kinetic equilibrium. 

This is quasi-equilibrium state since total number of particles is still approximately con- 
served, as can be seen from fig. [2] and triple interactions are not yet efficient. At the 
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FIG. 1: Depencence on time of energy densities of electrons (green), positrons (red), photons 
(black) and protons (blue) for initial conditions I. Total energy density is shown by dotted black 
line. Interaction between pairs and photons operates on very short timescales up to 10 -23 sec. 
Quasi-equilibrium state is established at t\ ~ 10 -14 sec which corresponds to kinetic equilibrium 
for pairs and photons. Protons start to interact with then as late as at ith — 10 -13 sec. 

moment t\ = 4 x 10~ 14 sec, shown by the vertical line on the left in fig. [3] and fig. HJ the 
temperature of photons and pairs is — 1.5, while the chemical potentials of these particles 
are ~ —7. Concentration of protons is so small that their energy density is not affected 
by the presence of other components; also proton-proton collisions are inefficient. In other 
words, protons do not interact yet and their spectra are not yet of equilibrium form, see 
fig. EJ The temperature of protons start to change only at 10~ 13 sec, when proton-electron 
Coulomb scattering becomes efficient. 

As can be seen from fig. HI the chemical potentials of electrons, positrons and photons 
evolved by that time due to triple interactions. Since chemical potentials of electrons, 
positrons and photons were negative, the particles were in deficit with respect to the thermal 
state. This caused the total number of these particles to increase and consequently the 
temperature to decrease. The chemical potential of photons reaches zero at t-i = 10~ 12 sec, 
shown by the vertical line on the right in fig. [3] and fig. HI which means that electrons, 
positrons and photons are now in thermal equilibrium. However, protons are not yet in 
equilibrium with other particle since their spectra are not thermal, as shown in the lower 
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FIG. 2: Depencence on time of concentrations of electrons (green), positrons (red), photons (black) 
and protons (blue) for initial conditions I. Total number density is shown by dotted black line. In 
this case kinetic equilibrium between electrons, positrons and photons is reached at ~ 10 -14 sec. 
Protons join thermal equilibrium with other particles at ith — 4 X 10 -12 sec. 

part of fig. [5j 

Finally, the proton component thermalize with other particles at 4 x 10 -12 sec, and from 
that moment plasma is characterized by unique temperature, 8 t h — 0.48 as fig. [3] clearly 
shows. Protons have final chemical potential v p ~ —12.8. 

This state is characterized by thermal distribution of all particles as can be seen from 
fig. [61 There initial flat as well as final spectral densities are shown together with fits of 
particles spectra with the values of the common temperature and the corresponding chemical 
potentials in thermal equilibrium. 

B. Case II 

We take the following initial conditions: power law spectral densities Eife) for protons, 
electrons and positrons with initial energy densities p p = 2.8 x 10 22 erg/cm 3 , p_ = 1.5 x 10 24 
erg/cm 3 , p + = 1.5 x 10 21 erg/cm 3 , respectively. We chosen flat spectral density for photons 
with p 7 = 2.8 x 10 24 erg/cm 3 . Initial baryonic loading parameter is set to B = 608, 
corresponding to a matter-dominated plasma, unlike the previous case. As in the case I, 
the most rapid reaction is electron-positron pair creation which starts to change the energy 
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FIG. 3: Depencence on time of dimensionless temperature of electrons (green), positrons (red), 
photons (black) and protons (blue) for initial conditions I. The temperature for pairs and photons 
acquires physical meaning only in kinetic equilibrium at t\ ~ 10 -14 sec. Protons are cooled by the 
pair-photon plasma and acquire common temperature with it as late as at tth — 4 x 10~ 12 sec. 

density of positrons at 10 -20 sec, see fig. [7J Initially most energy is in photons, followed by 
electrons and protons. In the course of the evolution the energy gets redistributed in such a 
way that in the final state most energy is transferred first to the electrons, then follow the 
protons, the photons and finally the positrons. In fig. [8] one can see that number densities 
of electrons and protons are almost equal with chosen heavy proton loading. Concentrations 
of particles almost do not change during evolution towards thermal equilibrium. 

Temperatures and chemical potentials of particles are shown in fig. 191 and [TU1 respectively 
Kinetic equilibrium is established at around 8 x 10~ 15 sec, marked by the vertical line. The 
temperature of pairs and photons at that moment is 9^ ~ 0.53, while the chemical potentials 
of these particles are V- ~ 1, i/ + ~ —0.9, z/ 7 ~ 0.1. Notice that chemical potentials of 
electrons and positrons are almost equal in magnitude and opposite in kinetic equilibrium, 
see fig. [TUl At this moment protons are not yet in equilibrium with the rest of plasma but 
already established kinetic equilibrium with themselves with the temperature 8 P ~ 0.18 and 
the chemical potential v v ~ —2. The common temperature is reached at the moment 10~ 13 
sec, which corresponds to thermal equilibrium. Final values of temperature is 8 t h — 0.47, 
while chemical potentials are v± ~ =pl, v v ~ —4.7. 
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FIG. 4: Depencence on time of dimensionless chemical potential of electrons (green), positrons 
(red), photons (black) and protons (blue) for initial conditions I. The chemical potential for pairs 
and photons acquires physical meaning only in kinetic equilibrium at ~ 10~ 14 sec, while for 
protons this happens at i t h — 4 x 1CP 12 sec. At this time chemical potential of photons has evolved 
to zero and thermal equilibrium has been already reached. 

The share of the proton energy density in the total energy density increased in course 
of time, see fig. 0, causing an increase in the baryonic loading parameter which reached in 
thermal equilibrium the value B = 780. 

Since concentration of protons is chosen to be large, proton-proton collisions become 
more important than proton-electron/positron collisions, in contrast to the case I. In fact, 
protons reached equilibrium temperature already at 10~ 16 sec, while they start to interact 
with electrons and positrons only at 1CT 15 sec. 

C. Case III 

We take the following initial conditions: the initial ratio between concentrations of elec- 
trons and protons is ? = n p /n^ = 10~ 3 . The total energy density is chosen in such a way 
that the final temperature in thermal equilibrium be 6th = 2. We set up flat initial spectrum 
for photons i£ 7 (e$) = const, and power law spectra for the pairs E±(e±) oc [e± — mc 2 ] 2 
and protons E p {e p ) oc [e p — Mc 2 ] 4 . Finally, the ratio of initial and final concentrations of 
positrons is chosen to be n + = I0~ 1 n t ^ 1 . Given these initial conditions the baryon loading 
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FIG. 5: Spectral density as function of particle energy for electrons (green), positrons (red), photons 
(black) and protons (blue) for initial conditions I at intermediate time moments t\ = 4 x 10 -14 sec 
(upper figure) and £2 = 10 -12 sec (lower figure). Fits of the spectra with chemical potentials and 
temperatures corresponding to thermal equilibrium state are also shown by yellow (electrons and 
positrons), grey (photons) and light blue (protons) thick lines. The upper figure shows the spectra 
when kinetic equilibrium is established for the first time between electrons, positrons and photons 
while the lower figure shows the spectra at thermal equilibrium between these particles. On both 
figures protons are not yet in equilibrium neither with themselves nor with other particles. 

parameter is B = 0.2. The initial conditions are chosen in a way to get larger temper- 
ature in thermal equilibrium, than in previous cases. Unlike the case II, the spectrum of 
protons is chosen steeper than the spectrum of pairs in order to make them colder in kinetic 
equilibrium. 

Equipartition of energies between pairs and photons occurs earlier than in the case I, at 
around 10~ 17 sec, see fig. [T2J since now concentrations of particles are higher. Concentrations 
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FIG. 6: Spectral density as function of particle energy are shown as before at initial and final 
moments of the computations. The final photon spectrum is black body one. 




FIG. 7: Depencence on time of energy densities for initial conditions II. Colors are as in the case 
I. Protons start to interact with other particles as late as at t ~ 10~ 16 sec. 

of pairs and photons equalize at 3 x 10~ 17 sec, see fig. [131 As in the case I, from this moment 
temperatures and chemical potentials of electrons, positrons and photons tend to be equal, 
see fig. [[Hand fig. [15] respectively, leading to kinetic equilibrium at around ~ 1CT 16 sec. 

At the moment t^, shown by the vertical line on the left in fig. [13] and fig. [15], the 
temperature of photons and pairs is #k — 2.2, the chemical potentials of these particles are 
i\ ~ —1.1, while the temperature of protons, having well established spectrum by this time, 
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FIG. 8: Depencence on time of concentrations for initial conditions II. Colors are as in the case I. 
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FIG. 9: Depencence on time of dimensionless temperature for initial conditions II. Colors are as 
in the case I. The pair-photon plasma is heating protons. Protons join thermal equilibrium at 
ith — 1CT 13 sec. 

is just 9 p ~ 0.09. 

Thermal equilibrium is reached in the electron-positron-photon plasma at around t t h — 
4 x 10~ 15 sec, shown by the vertical line on the right of fig. [HI and fig. [T3 Only at 4 x 10~ 14 
sec the temperature becomes common also with protons which are heated up during this 
time. The temperature at this final stage is 9 t h — 2 while the chemical potential of protons 
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FIG. 10: Depencence on time of dimensionless chemical potential for initial conditions II. Colors 
are as in the case I. The chemical potential of photons is almost zero in kinetic equilibrium. The 
chemical potentials of electrons and positrons are almost equal and opposite in kinetic equilibrium, 
to maintain electric neutrality. 
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FIG. 11: Initial and final spectral density as function of particle energy for initial conditions II. 
Fits of the final spectra with chemical potentials and temperatures are also shown. 

is v p ~ —33. Initial as well as final spectra are shown in fig. [T6j 

Since chemical potentials and temperatures approach their values in thermal equilibrium 
exponentially, i.e. ~ exp(— t/r c h. eq ), we determined the relaxation time constant r c h. eq for 
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FIG. 12: Depencence on time of energy densities for initial conditions III. Colors are as in the case 
I. Protons start to interact with other particles at about 10~ sec. 
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FIG. 13: Depencence on time of concentrations for initial conditions III. Colors are as in the case 
I. 



each of the cases considered from 

(24) 

where the fugacity for a given sort of particle is given by (|T3|) . Our results are presented in 
Tab. ED 
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FIG. 14: Depencence on time of dimensionless temperature for initial conditions III. Colors are as 
in the case I. Pairs and photons acquire the temperature at ~ 10 -16 sec. 




FIG. 15: Depencence on time of dimensionless chemical potential for initial conditions III. Colors 
are as in the case I. The chemical potentials equalize at ~ 10 -16 sec. 



VII. DISCUSSION AND CONCLUSIONS 



Results presented above clearly show the existence of two types of equilibrium: the kinetic 



and the thermal ones. Kinetic equilibrium in pair-photon plasma occurs when Ehlers 46 ] 
balance conditions (1B6I) . (1B9I) and ( 1B12I) are satisfied so that pair-creation, Compton and 
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FIG. 16: Initial and final spectral density as function of particle energy for initial conditions III. 
The spectrum of protons is chosen to be steeper than the one of electrons and positrons. Fits of 
the final spectra with chemical potentials and temperatures are also shown. 
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T ch'.eq» SeC 


2.2 x 1(T 13 


1.8 x 1(T 14 


9.5 x 1(T 16 


T c P h.eq> SeC 


6 x i(r 13 


1.8 x 1(T 14 


5.5 x 1(T 15 



TABLE III: Relaxation time constant for cases I-III. 

Bhabha/M0ller scattering processes all come in detailed balance. The electron-positron- 
photon plasma then is described by common temperature and nonzero chemical potentials 
which are given by (lB23l) .( lB24l) and ( 1B25I) . (1B26I) . Protons at this stage may or may not have 
yet established equilibrium with the spectrum ( fT5]) . depending on the value of the baryon 
loading parameter B. When B is small, as in the case I, proton-proton collision are inefficient 
since the rate (1D15|) is much smaller than (1D11I) . and the proton spectrum is shaped up by 
the proton-electron/positron collisions, reaching equilibrium form at a timescale given by 
(1D11I) . when other particles are already in thermal equilibrium. When B is large, as in the 
case II, protons have established their equilibrium temperature at a timescale given by ( 1D15I) . 
prior to the moment when kinetic equilibrium in the pair-photon plasma is established. 

As we have seen, the final spectra are completely insensitive to the initial spectra, chosen 
to be flat as in the case I, power-law as in the case II, or thermal ones. 
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The meaning of non-zero chemical potentials in kinetic equilibrium can be understood as 
follows. The existence of a non-null chemical potential for photons indicates the departure of 
the distribution function from the one corresponding to the thermal equilibrium. Negative 
value of the chemical potential generates an increase of the number of particles in order 
to approach the one corresponding to the thermal equilibrium state. Positive value of the 
chemical potential leads to the opposite effect, decreasing the number of particles. Then, 
since the total number of particles increases (or decreases), the energy is shared between 
larger (or smaller) number of particles and the temperature decreases (or increases). Clearly, 
as thermal equilibrium is approached, the chemical potential of photons tends to zero, while 
the chemical potentials of electrons and positrons are given by (1B27[) . to guarantee an overall 
charge neutrality. 

One of the basic assumptions in this work is that triple interactions are slower than binary 
ones, allowing to use reaction rates for triple interactions in kinetic equilibrium, explicitly 
depending solely on temperature, chemical potentials and concentrations of particles. For 
pure electron-positron plasma in the range of energies of interest (CEJ) there is a hierarchy 
of relevant timescales: binary interactions are clearly faster than triple ones. However, 
when protons are also present, the proton-proton timescale may be shorter or longer than 
the corresponding binary interactions timescales for the pure pair plasma. This violates 
our assumption and therefore leads to loss of quantitative accuracy, although still keeping 
qualitative results valid. In order to overcome this difficulty and produce quantitatively 
precise results exact QED matrix elements must be used for calculation of emission and 
absorption coefficients. 

Notice that there is some discrepancy between our final spectra and their thermal fits for 
high energy. This is due to poor energy resolution with adopted grid. The result converges 
with higher resolutions, but it is limited by the available computer memory. In addition, 
the code is quite time-consuming and processor time increases with number of operation as 
third power of the number of energy intervals. 

In order to resolve proton-electron/positron scattering the number of energy intervals 
should be increased as M/m comparing to the case of pure pair plasma. Even using in- 
homogeneous energy grid with uniform energy step up to the peak of the spectrum dp/de 
and decreasing energy step as e~ x for higher energies, we have obtained acceptable results 
with about 10 3 intervals for this reaction. Using such fine grid is impossible in practice. 
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On the other hand, a small parameter m/M expansion can be adopted. In this way we 
have introduced the mass scaling, described in Appendix [Gl which gives quite good accu- 
racy for about 10 2 intervals in energy with inhomogeneous grid, described above. Finally, 
it is important to stress that our code allows solution of the Boltzmann equations for long 
time intervals and timescales, which may differ up to 10 orders of magnitude, from electron- 
positron creation and annihilation process up to proton-electron/positron scattering, see fig. 



El unlike approaches based on Monte-Carlo technique [45|] . This gives us the possibility to 
follow thermalization process up to reaching steady solution, i.e. thermal equilibrium. 

The assumption of the constancy of the energy density is only valid if the following three 
conditions are satisfied: 

• plasma is optically thick for photons. This leads to the constraint on the spatial 
dimensions R ^> (n th cr T ) 1 ~ 1CT 5 cm. 

• neutrinos are not produced. This gives the constraint on the temperature 9 <C 7 x 10 2 . 

• plasma does not expand. Given t dyn = (j^^f ) 1 3> tth, this leads to R ^> 1CT 2 cm. 

To summarize, we have considered the evolution of initially nonequilibrium optically thick 
electron-positron-photon plasma with proton loading up to reaching thermal equilibrium on 
a timescale tth < 10 -11 sec. Starting from arbitrary initial conditions we obtain kinetic 
equilibrium, on a timescale < 10~ 14 sec, from first principles, solving numerically the 
relativistic Boltzmann equation with collisional integrals computed from exact QED matrix 
elements. 

The general theoretical framework here presented can be further applied by consider- 
ing thermalization of different relativistic particles predicted by extensions of the standard 
model of particle physics with the lepton plasma in the early Universe. The occurence of 
thermalization process of electron-positron plasma in GRBs on a much shorter timescale 
than the characteristic acceleration time [48] is crucial. Such acceleration timescales are 
indeed sharply bounded (shorter than 10 3 sec in the laboratory frame). Determination of 
thermalization timescales as functions of the relevant parameters is important for high en- 



ergy plasma physics 49(, 50j. Finally, these results can in principle be tested in laboratory 



experiments aiming the generation of electron-positron pairs. 



26 



We thank the anonymous referee for comments which improved the presentation of our 
results. 



APPENDIX A: CONSERVATION LAWS 

Conservation laws consist of baryon number, charge and energy conservations. In addi- 
tion, in binary reactions particle number is conserved. 

Energy conservation law can be rewritten for the spectral density 

sl> = °. or |En. = o, (ad 

i 

where 



/■ei ]W +Aei ]W /2 

Y^= / Eide. (A2) 



Particle's conservation law in binary reactions gives 

!5>=°. « |E£ = o- (A3) 

Since baryonic number is conserved, therefore the number density of protons is a constant 

dn r 



2 = 0. (A4a) 



For electrically neutral plasma considered in this paper charge conservation implies (TTT1) . 

APPENDIX B: DETERMINATION OF TEMPERATURE AND CHEMICAL PO- 
TENTIALS IN KINETIC EQUILIBRIUM 

Consider distribution functions for photons and pairs in the most general form (IT51) . If 
one supposes that reaction rate for the Bhabha scattering vanishes, i.e. there is equilibrium 
with respect to reaction 

e + + e~ <-> +e +/ + e~', (B5) 
then the corresponding condition can be written in the following way 

f + (l - /+')/_(! - /_') = /+'(! - /+)/_'(! + /_), (B6) 
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where Bose-Einstein enhancement along with Pauli blocking factors are taken into account 
for generality, it can be shown that electrons and positrons have the same temperature 



6+ = 6- = 6±, (B7) 



and they have arbitrary chemical potentials. 

With (IB 70 analogous consideration for the Compton scattering 

e ± + 7 <-> +e ±/ + 7', (B8) 

gives 

/±(i - /± ; )A(i + A) = /±(i - /±)A'(! + A)> ( B9 ) 

and leads to equality of temperatures of pairs and photons 

0± = # 7 = 9 k , (BIO) 

with arbitrary chemical potentials. If, in addition, reaction rate in the pair-creation and 
annihilation process 

e ± + e T 4-^7 + 7' (Bll) 

vanishes too, i.e. there is equilibrium with respect to pair production and annihilation, with 
the corresponding condition, 

+ / 7 )(1 + A) = AA(i - A)(i - /-), (B12) 

it turns out that also chemical potentials of pairs and photons satisfy the following condition 

v+ + v- = 2u 7 . (B13) 



However, since, generally speaking, v n 7^ the condition (1B13|) does not imp 



y 



51] , and we 



These considerations were for the first time applied by Ehlers in [46], see also 
will call flB6j) . flB9l) and (IB12I) the Ehlers balance conditions. 

Analogous consideration for the detailed balance conditions in different reactions lead to 
relations between temperatures and chemical potentials summarized in table IIVI 

The timescales of pair production and annihilation processes as well as Compton scat- 
tering are nearly equal in the range of energies of interest and are given by (TL4"|) . Therefore, 
kinetic equilibrium is first established simultaneously for electrons, positrons and photons. 
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Interaction 


Parameters of DFs 


I 


e + e~ scattering 


9 + = 9-, Vv + ,v_ 


II 


e ± p scattering 


9 P = 9±, Vv±,v p 


III 


e ± 7 scattering 


9*y = 9±, \/f~J,U± 


IV 


pair production 


v+ + = 2i>y , if 6ty = 9± 


V 


Tripe interactions 


Ury, v± = 0, if 9y = 9± 



TABLE IV: Relations between parameters of equilibrium DFs fulfilling detailed balance conditions 
for the reactions shown in Tab. HI 



They reach the same temperature, but with chemical potentials different from zero. Later 
on, the temperatures of this electron-positron-photon plasma and the one of protons reach 
a common value. 

In order to find temperatures and chemical potentials we have to implement the following 
constraints: energy conservation (lAlj) . particle number conservation (1A3j) . charge conserva- 
tion ffTTj) . condition for the chemical potentials (IB13I) . 

Given (|T5l) we have for photons 

1 

7 



th 



36. 



(B14) 



for pairs 



and for protons 



P± 



n±mc* 



Pv 



Mn p c 2 



i + --e p , 



1 /WM\ 3/2 

^^VHw exp 



u p - 



M 



Or, 



e 



p ■ 



(B15) 

(B16) 
(B17) 



where we assumed that protons are nonrelativistic, and we denoted the Compton volume 

3 







2nh 



7T \ mc 



(B18) 
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and functions ji and j 2 are defined as 
j 1 {6) = eK 2 {6- 1 )- 



fe-^ 3 / 2 , 0^0 



20 3 



30 



h{8) 



2 . 0^0 

30, ^ oo 



(B19) 



(B20) 



For pure electron-positron-photon plasma in kinetic equilibrium, summing up energy 

densities in flBTij) . dBTB]) and using (iBTD. dBTOj) and (1BT31) we obtain 

2mc 2 ( z/k 
Ok. 



e+,e — ,7 



Vn 



exp 0£\ [Mt+hie*)]^ 



(B21) 



and analogously for number densities we get 



E * = £«i>(£)[«i+AW. 



(B22) 



From flB21j) and (IB22[) two unknowns, and 0k can be found. 

When protons are present, in most cases the electron-positron-photon plasma reaches 
kinetic equilibrium first, while protons join the plasma later. In that case, the temperature 
of protons 9 P is different from the rest of particles, so while 9 + = 0_ = 7 = 0k, 9 P ^ 0k- 

Then summing up energy densities in (1B14j) . (1B15l) we obtain 



e+,e",7 



mc 



x 60k exp 



Et 
k 



_ n p Vo 
Ji(0k) 

2j 1 (0 fc )exp 



exp 



E± 
k 



and analogously for number densities we get 



E 



rij 



e T ,e ,7 

x 60k exp 



1 

V 



n p V 



exp 



E± 
k 



+ 2j 1 (0 k )exp 



k 

n p V 



"0 k 

E+ 
k 



x 



(B23) 



(B24) 



From ( 1B23I) and (1B24I) two unknowns, v + and 0k can be found. Then the rest of the 
chemical potentials are obtained from 



exp 



0k" 



exp 



0,< 



n p Vo 
Ji(0k) 



<>X1 * 1 4 



■ u + 

ex P 7T 



n p V 
Ji(0k) 



ex P "TT 



(B25) 
(B26) 
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The temperature and chemical potential of protons can be found separately from 
dBHD,dB17D. 

In thermal equilibrium z/ 7 vanishes and one has 

u+ = -u_, (B27) 

which both reduce to z/_ = v + = for n p = 0. At the same time, for n p > one always has 
V- > and v + < in thermal equilibrium. The chemical potential of protons in thermal 
equilibrium is determined from ( 1B17I) for 9^ = 6^, where 6 t h is the temperature in thermal 
equilibrium. 



arcsinh 



n p Vo 



2ji 



APPENDIX C: BINARY INTERACTIONS 



1. Compton scattering 7e ± — ► j'e^' 

The time evolution of the distribution functions of photons and pair particles due to 
Compton scattering may be described by [341]. [52 1 



dt 

df±(p,t) 



dt 



= /dk'dpd P V Wk ^ ;k , p [/ 7 (k',t)/ ± (p',t)-/ 7 (k,t)/ ± (p,t)], (CI) 

y'p±' J 

dkdk'dp>Vw k ,, p , M Mk>,t)f ± (p',t) - f 7 (k,t)f±(p,t)}, (C2) 



where 



2^6 



»k',p';k,p 



(2n) 2 V 

is the probability of the process 



5(e 7 - e± - e' - e' ± )S(k + p - k' - p' 



|M /4 | 2 



16e 7 e- t e 7 e / ± ' 



(C3) 



|M 



fi\ 



2\ 2 a 2 



m 2 c 2 
s — m 2 c 2 



+ 



m 2 c 2 
u — m 2 c 2 



+ 



m 2 c 2 
— m 2 c 2 



+ 



m 2 c 2 
u — m 2 c 2 



m 2 c 2 



a 



m 2 c 2 



(C4) 



4 \ u — m 2 c 2 s — m 2 c 2 _ 

is the square of the matrix element, s = (p + 6) 2 and u = (p — t') 2 are invariants, t = 
(e 7 /c)(l,e 7 ) and p = (e±/c)(l, /3±e±) are energy-momentum four vectors of photons and 
electrons, respectively, dp = de±doe 2 t /3±/c 3 , dk' = de'^e' 2 do'^ / (? and do = dfidcj). 
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The energies of photon and positron (electron) after the scattering are 
, _ e±e 7 (l-/4b±-b 7 ) _ 



b t = Pi /p, b\ = p'/p>, b' ± = (/3±e ± b± + e 7 b 7 - e;b;)/(^4). 

For photons, the absorption coefficient (flOl) in the Boltzmann equations (HI) is 



r*-W*»t Udf,Y hs f. ,, T ^\M fl \ 2 h 2 c 2 f 
Xl 7 fj = — I J = / dn±do J cs ^— ; — / 7 , (C6) 



c V ^ / 7e ±_y e ±' </ 7 16e±e 7 e' ± 

where <inj = deidoiefPifi/c 2, = dejdoj£/j/(27rej). 

From equations ( ICll) and ( 1C6I) , we can write the absorption coefficient for photon energy 
density E n averaged over the e, /x-grid with zone numbers to and k as 



7 e±_» 7 ' e ±' 



/7 

-7,0; </<E 7 eAe 7 

1 /" , dlM/ilW 
= -7 / dn y dn±do Jcs-!— , (C7) 

where the Jacobian of the transformation is 

e' e' 

e 7 e± (1 -/7±b 7 -b ± ) 

Similar integrations can be performed for the other terms of equations (IClfl . (1C2I) . and we 
obtain 



e ,2 \M f A 2 h 2 c 2 



<^ 7 ' e± ' = T / d ni dn ± do\J cs 7 Jil •; : , (C9) 

Ae 7jt4 , J^eA^, u 7 16e±e 7 e^ 



e ±^ ' ±' 1 /" / e' \M fi \ 2 h 2 c 2 
rf± u) ^ e =- / dn 1 dn±do J cs - L — , (CIO) 



1 y . . . , . e 7 |M /4 |W 



(xE)±^' e = -z / dn 1 dn±do'J cs - :L — — . (Cll) 

Ae ±iW 7 e±eA e±, w 16e 7 e ± 

In order to perform integrals (10T[) - (101 1|) numerically over <fi (0 < <fi < 2n) we introduce a 
uniform grid <^ T i/ 2 with 1 < I < Z max and Afa = (0/+1/2 - <j>i-i/2)/2 = 2n/l max . We assume 
that any function of in equations (1C7l) - (lC9j) in the interval A(pj is equal to its value at 
4> = 4>j = (</>/-i/2 + </>z+i/2)/2. It is necessary to integrate over only once at the beginning 
of calculations. The number of intervals of the 0-grid depends on the average energy of 
particles and is typically taken as / max = 2fc max = 64. 
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2. Pair creation and annihilation 7172 ^ e e + 

The rates of change of the distribution function due to pair creation and annihilation are 



df l3 {k h t) 
dt 



7172— > e _ e _t 



/ 



dk j dp-dp + Vw p _ i p + .u 1 ,k 2 f'riO s -i,t)fr r2 ('k-2,t) , (C12) 



\ / e~ e+— >7i72 J 



(C13) 



e _ ef — >7i72 

for i = 1, j = 2, and for j = 1, i = 2. 



' df ± ( P± ,t) 
dt 



J = / dp T dk 1 dk 2 Vw p _ tP+ - ]lu]l J 1 (k ll t)f y (k 2 ,t) , (C14) 

/ ^^2— >e — e+ ^ 



( ^^'^ J = ~ [ dp T dk 1 dk 2 Vw klM2 . p _ jP+ f^{p^,t)f + (p + ,t) , (C15) 

\ °^ / e-e+^7172 ^ 



where 



^..p+jkLka = x 2T/ ^( e - +£+-£!- e 2 )5(p_ + p+ - ki - k 2 ) — . (C16) 

Here, the matrix element |M/j| 2 is given by equation (1C4I) with the new invariants s = 



(p_ - h) 2 and u = (p_ - t 2 ) 2 , see |41| 



The energies of photons created via annihilation of a e pair are 

m 2 c 4 + e_e + (l - /3_/?+b_-b+) 

61 bl = ~T\ a u u vT — n a k k V ' e 2 (bi) = e_ + e+ - ei , (C17) 

e_(l - /L.b_-bi) + e+(l - p+b+'bi) 

while the energies of pair particles created by two photons are found from 



. B =F VB 2 - AC . 
e-(b-) = j , e + (b_) = ei + e 2 -e_, (C18) 

where A = (e l + e 2 ) 2 - [(eibi + e 2 b 2 )-b_] 2 , B = (e l + e 2 )eie 2 (l - b r b 2 ), C = m^Keibj + 
e 2 b 2 )-b_] 2 + e 2 e2(l — bi-b 2 ) 2 . Only one root in equation ( 1C18I) has to be chosen. From 
energy-momentum conservation 

t 1 + t 2 -p_=p + , (C19) 
taking square from the energy part we have 

e 2 + e 2 + e 2 + 2 ei e 2 - 2 ei e_ - 2e 2 e_ = e 2 (C20) 
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and taking square from the momentum part we get 

4 + 4 + <?-P- + 2eie 2 b r b2 - 2e 1 e_0_b r b_ - 2e 2 e_/3_b 2 -b_ = (e+/3+) 2 . (C21) 

There are no additional roots because of the arbitrary e + 

eie 2 (l - b r b 2 ) - eie_(l - /3_b r b_) - e 2 e_(l - /3b 2 -b_) = 0, (C22) 
e_/?„(eibi + e 2 b 2 )-b_ = e_(ei + e 2 ) - eie 2 (l - b r b 2 ). 

Eliminating /3 we obtain 

e 2 e 2 (l - b r b 2 ) 2 - 2 £l e 2 (l - b r b 2 )( ei + e 2 )e_ + {( £l + e 2 ) 2 - [(e^ + e 2 b 2 )-b_] 2 } e 2 _ = 
= [( ei b 1 + e 2 b 2 )-b_](-m 2 ), (C23) 

Therefore, the condition to be checked reads 

e_/3_ [(exbi + e 2 b 2 )-b_] 2 = [e_( ei + e 2 ) - ( ei e 2 )(l - b r b 2 )] [( Cl bi + e 2 b 2 )-b_] > 0. (C24) 

Finally, integration of equations (1012[) - (1C15[) yields 

16e_e + e 2 Je 2 eAe^, u 16e_e+ 

(C25) 



(x*0£f 7172 = a — / ^V4 '' J " + / rfV J ( 



Ae e>£J yJe-eAe^ 16e+e 2 Je+eAe e , u 16e_e 2 

(C26) 



(C27) 



^e e ,o; We_eAe e ,o, 106ie 2 e + Je+&Ae e ^ 10€ie 2 

(C28) 

where d 2 n± = dn-dn + do\ 1 d 2 n 1 = dn^dn^do-, dn± = de±do±e 2 ±(3±f± 1 dn 112 = 
dei^doi^4 2/71,2 an d the Jacobian is 

(e+ + e_) /3_ - (eibi + e 2 b 2 ) -b_ 
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3. M0ller scattering of electrons and positrons efef — > efef' 



The time evolution of the distribution functions of electrons (or positrons) is described 



by 



( dfi{pj,t) 
V dt 



dpjdpidp' 2 Vw 



Pi:P' 2 ;Pl,P2 



[/l(p , l,t)/ 2 (p / 2 ,t)-/l(Pl,t)/2(p 2 ,t)], 

(C30) 



with i = 1, j = 2, and with j = 1, z = 2, and where 
£ 2 c 6 



w PijPa;pi>Pa 



(2tt) 2 ^ 



+ e 2 - e; - e 2 )5( P i + p 2 - pi - p 2 ) Wfl ^ 



16eie 2 e / 1 e / 2 ' 



\M fi \ 2 = 2 V 



1 

+— 

u 



t 2 



r s 2 + t 2 



S 2 + U 2 2 2/. 2 2\ 

h 4m c (t — m c ) 



+ 



+ 4m c (u — m c 



— f m 2 c 2 ) ( 3m z c 



tu V2 



,2 2 



(C31) 
(C32) 
(C33) 



with s = (pi + po} 2 = 2(m 2 c 2 + pip 2 ), t=(p l - p[) 2 = 2(m 2 c 2 - pip'J, and u=(p 1 -p 



I \2 



2(m 2 c 2 -pip' 2 ) [41]. 

The energies of final-state particles are given by ( 1C18I) with new coefficients A = 
(ei + e 2 ) 2 - (ei/3ibrb; + e 2 /3 2 b 2 -b' 1 f , B = (e 1 + e 2 )[m 2 c 4 + e x e 2 (l - fr&bxba)], and 
C = m 2 c 4 (e 1 /3 1 b r b , 1 + e 2 f3 2 b 2 -b' 1 ) 2 + [m 2 c 4 + ei e 2 (l - /3 1( 5 2 b r b 2 )] 2 . The condition to be 
checked is 

+ e 2 ) - m 2 c 4 - ( ei e 2 )(l - /3i/3 2 b r b 2 )] [(eiftbi + e^b,)-^] > 0. (C34) 



Integration of equations (IC30j) . similar to the case of Compton scattering in Section |C II 
yields 



eie 2 — *e' 1 e' 2 
le,ui 



Ae f 



d nJ n 



(XE)1 



Ae 



d nJ„ 



s' 2 P[\M fl \ 2 h 2 ^ 
16eie 2 e' 2 



+ 



16eie 2 



(C35) 



+ / d nJ ms - 

16e 2 e 2 y e2 6Ae e , u 16eie 2 



where <i 2 n = dnidn 2 do[, dni )2 = de\^ 2 doi )2 e\ 2 j3i )2 f 121 and the Jacobian is 

e' 2 {3' 2 



(e' 1 + e' 2 )P' 1 -(e 1 p 1 b 1 + e 2 p 2 b 2 )-b , 1 



(C36) 



(C37) 
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4. Bhaba scattering of electrons on positrons e e + — > e e 

The time evolution of the distribution functions of electrons and positrons due to Bhaba 
scattering is described by 

f df ± (p ± ,t) \ = r d ^ dpj >_ dp , +VWpi , pi p _ iP+ [/_(p'_,t)/ + (p' + ,t)-/_(p_,t)/ + ( p+ ,t)], 

V 01 J e -e+^e-'e+' J 

(C38) 

where 

w P '_,p; ; p_,p + = ( 2n) 2 V 5 ( e - + e + ~ e - " e '+)<KP- + P+ - P- - P+) 16e _ e \>_ e > ' ( C39 ) 



and |M/i| is given by the equation (1C33I) . but the invariants are s = (p_ — P+) 2 , t = (p+— P_ 



+ J 2 



'unctions of the outgoing particle directions 



4l| 



and u = (p_ + p + ) 2 . The final energies e'_, e' + are 
in a way similar to that in Section IC 3\ see also 
Integration of equations (1C38I) yields 

V±£ = T / dn±Jbs — 77 } + / d 2 n ± J bs — ^ 

(C40) 

rv ^ e - e+ _-' e+ ' _ 1 ( [ », j ^>-\M fl \ 2 h 2 c 2 f e>MMgW 
KX&)±,u - ~r / d n±As — ; h / a n ± J hs — ; 

(C41) 

where d 2 n' ± = dn_dn + do'_, dn± = d€±do±e 2 t /3±f±, and the Jacobian is 

e[ 13' 

Jhs = (Z + e' + )(3'_ - (eXh- + e + (3 + b + )-b'_ ' ( ° 42) 



Analogously to the case of pair creation and annihilation in Section ( IC 21) the en- 



ergies of final state particles are given by (1C18I) with the coefficients A = (e_ + 
e+) 2 - (e_/3_b_-b'_ + e+/3 + b + -b'_) 2 , B = (e_ + e+) [m 2 c 4 + e_e + (l - /3_/3+b_-b + )], C = 
[m 2 c 4 + e_e+(l - / 3_/3+b_-b+)] 2 + m 2 c 4 [e_/3_b_-b'_ + e + /3+b + -b'_] 2 . In order to select the 
correct root one has to check the condition (1C34j) changing the subscripts 1 — ^ — , 2 — > -h. 
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APPENDIX D: BINARY REACTIONS WITH PROTONS 



1. Compton scattering on protons jp — ► j'p' 

The rate for this process t~£, compared to the rate of Compton scattering of electrons 
t~l is much longer, 

^ = l> -^- (D1 » 

Moreover, it is longer than any timescale for binary and triple reactions considered in this 
paper and thus we exclude this reaction from the computations. 



2. Electron-proton and positron-proton scattering e±p — > e'±p' 

The time evolution of the distribution functions of electrons due to ep — > e'p' is described 



by 

' df±(p,t) 
dt 

df P (<j,t) 
dt 

where 



dqdp , dqVw P '.q' ; p l q[/±(p / J t)/p(q / ,t) - /±(p,*)/ P (q,f)], (D2) 
dpdp , dqV«;p/, q ,.p, q [/ ± (p , ,t)/ p (q , ,t) - /±(p,*)/„(q,*)], (D3) 



h 2 c 6 „, , x \M fi 



2 



W *«™ = J^fV 5 ^ + e *~ < - + q - P' - q/) 16^' (D4) 

\M fl \ 2 = 2 6 7rV^ j^(s 2 + « 2 ) + (mV + M 2 c 2 )(2t - m 2 c 2 - M 2 c 2 )| , (D5) 

the invariants are s = (p + q) 2 = m 2 c 2 + M 2 c 2 + 2p • q, £ = (p — p') 2 = 2(m 2 c 2 — 
p • p') = 2(M 2 c 2 - q • q') and u = (p - q') 2 = m 2 c 2 + M 2 c 2 - 2p • q', s + t + u = 
2(m 2 c 2 + M 2 c 2 ). The energies of particles after interaction are given by flC18j) with 



A = (e± + e p ) 2 - [(e±/3±b ± + e p f3 p b p )-b' ± ] \ B = (e± + e p )[m 2 c 4 + e±e p (l - /3±/3 p b±-b p )], 
(7 = m 2 c 4 {(e±/3±b±-b± + ep/^bp-b^) 2 + [m 2 c 4 + e±e p (l — /3±/3 p b±-b p )] } 2 . The correct root 
is selected by the condition f]C34|) with the substitution 1 — > ±, 2 — > p. 
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Absorption and emission coefficients for this reaction are 



1 f ^ 7 



(X^)±,^ = T / dn±dn p do'±J ep — — , (D6) 



m ^ = U„ dnMJ - ■ (D7) 



i / ,„ ^ 7 4/44|M„iW 



^ = -X / dn±dn p do±J ep — — , (D8) 

1 f , e' 2 (3Le'\M fi \ 2 h 2 c 2 , x 

= ^~ / dn ± dn p do' ± J ep ^/ 7 ; , (D9) 

where drij = deidoieffiifi, i = ±,p, and the Jacobian is 

e' 

6P (4 + e' p W ± - (e p /3 p h p + e ± (3 ± b ± ).b' ± ' { ] 

The rate for proton-electron (proton-positron) scattering is 

e ± <e p . (Dll) 

3. Proton-proton scattering p\P2 — ► P1P2 

This reaction is similar to eie 2 — > e^e^, described in Section |C 31 The time evolution of 
the distribution functions of electrons is described by 

( dfi{pi,ty 



V dt 



PiP2-»piP2 



dqjd^d^Vw^.^^ [/i(qi, t)/ 2 (q£, t) - / x (qi, t)/ 2 (q 2 , *)], 

(D12) 

with j = 3 — i, and where 

^ S(e 1 + e 2 - e; - 4)%! + q 2 - ^ - q 2 ) J^y, (D13) 



Qi,q 2 ;qi,q 2 - ( 2vr )2^1 ' ^ ^1 ^VHl , Hi Ml H 2 , £ , , 



s 2 + u 2 



+ AM 2 c 2 (t- M 2 c 2 ' 



2 



+ 



u 



f t2 +4M 2 c 2 ( M -MV) 



+ ^(l- MV )(i- 3MV )^ ( di4 ) 



2 

and the invariants are s = (q 1 + q 2 ) 2 = 2(M 2 c 2 + qi ■ q 2 ), t = (q 1 — q' x ) 2 = 2(M 2 c 2 — qi ■ q[ 
and u = (qi - q^) 2 = 2(M 2 c 2 - qiq' 2 ). 
For the rate we have 



_ i / 171 Tin i / ?7l 

Vp-\/m-^, ^Vm v± ' v± " c - (D15) 
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APPENDIX E: THREE-BODY PROCESSES 



We adopt emission coefficients for triple interactions from 31]. 
Bremsstrahlung 



' 3 e \ mc 



2 \ 2 



hi 



4f(11.2 + 1O.40 2 



16 2ac ( e 



n>n_ 



3 e \ mc 



.2 \ 2 



In 



4f (1 + 1O.40 2 ) 







2 \ 2 



In 



4f (1 + 3.420) 



91 \\f29 + 20 2 
'ej exp(l/0)K 2 (l/0) ! 

\/2 + 20 + 29 2 
exp(l/9)K 2 (l/ey 

1 + 20 + 29 2 
exp(l/9)K 2 (l/9Y 



(El) 

(E2) 
(E3) 



' 3 e \mr 

where £ = e -0 - 5772 , and K2(l/9) is the modified Bessel function of the second kind of order 
2. 

Double Compton scattering 



2 \ 2 



2 



3 e UcV 1 + 13.910 + 11.O50 2 + 19.920 3 ' 



Three photon annihilation 



r^ ±e=F ^ 77 7 = n + n_ac 



2 \ 2 1 !(21n 2 2^ + ^-| 



£40 + ^(21n 2 2£0 + ^-±) 



where we have joined two limiting approximations given by 
Radiative pair production 



3l| 



(E4) 



(E5) 



77'— +7"e ± eT 



V, 



Electron-photon pair production 



e^e^- >77'7" 



IK 



n+ri- 



K 2 {l/9) 
29 2 



± ±1 + m 

h 



(n+ + n-)n 7 ac(^) exp (-|) 16.10 



30.541 



9 < 2, 



Proton-photon pair production 



P7— >p'e ± e = F 



y ac (i) exp (-§) TT i^, < 1.25277, 
n^ac (^) 2 [f (ln2£0 + 1.7) - f ] , > 1.25277. ' 

We use the absorption coefficient for three-body processes written as 



(E6) 



(E7) 



(E8) 



(E9) 
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where 7/ 3p is the sum of the emission coefficients of photons in the three particle processes, 
E° q = 27re 3 /^ q /c 3 , where f° q is given by (fl5j) . 

From equation (1221) . the law of energy conservation in the three-body processes is 

/ $>, 3p - x- P £,)^ = . (E10) 

For exact conservation of energy in these processes we introduce the following coefficients of 
emission and absorption for electrons: 



X 3 P = I(v 3 y p -X 3 , p E,)ded^ 
' f E e dedfi 

or 



e 3p = 0, J (rff - X^Ejdedfi > , (Ell) 



3 P j (r? 3 P _ x ^ Ei)ded ^ 



V, 



6 " ^ P = 0, J(v?-xfE j)d ed,<0. (E12) 



APPENDIX F: CUTOFF IN THE COULOMB SCATTERING 

Denote quantities in the center of mass (CM) frame with index 0, and with prime after 
interaction. Suppose we have two particles with masses mi and m 2 . The change of the angle 
of the first particle in CM system is 

6»io = arccos(bi -b' 10 ), (Fl) 

the numerical grid size is A# g , the minimal angle at the scattering is # m i n . 
By definition in the in CM frame 

Pio + P20 = 0, (F2) 

where 



(r-i)(NpO-r-^ 

c c 



N, i = l,2, (F3) 



PiO =Pi + 

and 

€i = r(e i0 + Vp«,). (F4) 
Then for the velocity of the CM frame we have 

X = C P1±^, N V r = _^. (F5) 
c ei + e 2 \ . I x _ 
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By definition 



and then 



'10 



>20, 



b' 10 = b'. 



20- 



(F6) 



IP10I 



P20I = Po 

1 



c 



where 



e 20 wi 2 c 4 , 



( 6l + e 2 ) 2 -r 2 (m 2 -m 2 )c 4 

2(e! + e 2 )r 
(6 1 + 6 2 ) 2 + r 2 (m|-mf)c 4 

2(ei + e 2 )r 



(F7) 



eio 



e 2 o 



(F8) 
(F9) 



Haug [44| gives the minimal scattering angle in the center of mass system 

"mi 



(F10) 



-McD ( 7r + l)y/2( lr -1) 

where Ai, as above, is the reduced mass, the maximum impact parameter (neglecting the 
effect of protons) is 

° 2 Po (Fll) 



D 



u eio 



and the invariant Lorentz factor of relative motion (e.g. [4J]) is 



1 



7r 



eie 2 - pip 2 <r 
mim 2 c 4 



(F12) 



In the CM frame we finally obtain 



mc 



(l - /9 2 cos 6» min ) 



Since it is invariant, we then replace t in the denominator of |M/j| 2 in (IC33I) by the 
value t^Jl + t^Jt 2 to implement the cutoff scheme. Also at the scattering of equivalent 
particles we remove the case of exchange of particles as well as scattering on small angles, 
in other words we change u in the denominator of |M/j| 2 in (1C33l) .( lD5l) and (1D14I) by the 



value u^/l + t^Ju 2 . 
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APPENDIX G: MASS SCALING FOR THE PROTON-ELECTRON/POSITRON 
REACTION 



Since proton mass is larger than electron mass-energy M 3> m, e then for the CM frame 

V^^A r«l, Jxwl, (Gl) 
e[ - ei w V (e^ - e 01 )p oc (G2) 

and also 



s 2 



c 4 

W 2 



C 



M 4 + 4mM 3 + 6m 2 M 2 , (G3) 

. - M 4 — 4mM 3 + 6m 2 M 2 , (G4) 

M /i | 2 oc4(6m 2 -2t)M 2 , (G5) 

6 



while 



t = ~ 2m ^ (1 - ee ° e - o) = (G6) 



2 

eO 



= ^fi^ ^ [l + QlM- 1 )] (G7) 

for small angles. 

This leads to the following scaling for the reaction rate 

We can therefore calculate t^ , {xE)f° for a pseudo-particle with mass M 3> m, e 
instead of M and obtain 



VZ « (G9) 

(x£): p . - ^ (x^)r • (gio) 

For such purpose we selected the mass of this pseudo-particle as M = 20m. 
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APPENDIX H: THE DEFINITION OF MATRIX ELEMENTS 



Following 41( define the scattering matrix, being composed of real and imaginary parts 

S fl = 5 fi + i (2irh) 4 5 (4) (p y - Tfi, (HI) 

where 5fi is the unity matrix, 5^ stands for the four-momentum conservation and the 
elements of Tfi are scattering amplitudes. 

The transition probability of a given process per unit time is then 



w fl = c(2Tihf5^ (p f -pi)\T fi \ 2 V, 



(H2) 



where V is the normalization volume. 

For a process involving a outgoing particles and b incoming particles the differential 
probability per unit time is defined as 



dw = c(2vr/i) 4 5 (4) (p f - pi) \M fi \ 2 Vx 



(H3) 



x 



n 



he 



n 



dp' a he 



where and e' a are respectively momenta and energies of outgoing particles, e& are energies 
of particles before interaction, are the corresponding matrix elements, 5^ stands for 
energy-momentum conservation, V is the normalization volume. The matrix elements are 
related to the scattering amplitudes by 



M 



fi 



n 



he 



n 



he 
2ZV 



Tfi- 



(H4) 



For a binary process with 2 incoming and 2 outgoing particles it is convenient to introduce 
the differential cross-section. In fact, the differential probability for incoming particles with 
four momenta pi and p 2 , energies t\ and e 2 and masses m 1 and m 2 respectively, is just the 
product of the differential cross-section and the flux density 



where 



J 
I 



dw = jda, 



cl 



eie 2 V 



mi17l2C z 



(H5) 

(H6) 
(H7) 
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2 

In the CM reference frame the relation between the cross section and \Mfi\ acquires simplest 



form if cross-section is independent on the azimuth of p[ relative to pi then 

h 2 c A . , , l2 dt . . 

^=64^7' (H8) 

t = (Pi -p 2 ) 2 , (H9) 

cft = 2|pi||p' 1 |dcostf, (HIO) 

where d is the angle between pi and p^. 
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